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Abstract 

The energy spectrum and primary composition of cosmic rays with energy between 
3 X 10^^ and 3 x 10^*^ eV have been studied using the CASA-BLANCA detector. 
CASA consisted of 957 surface scintillation stations; BLANCA consisted of 144 
angle-integrating Cherenkov light detectors located at the same site. CASA mea- 
sured the charged particle distribution of air showers, while BLANCA measured the 
lateral distribution of Cherenkov light. The data are interpreted using the predic- 
tions of the CORSIKA air shower simulation coupled with four different hadronic 
interaction codes. 

The differential flux of cosmic rays measured by BLANCA exhibits a knee in the 
range of 2-3 PeV with a width of approximately 0.5 decades in primary energy. The 
power law indices of the differential flux below and above the knee are —2.72 it 0.02 
and —2.95 ± 0.02, respectively. 

We present our data both as a mean depth of shower maximum and as a mean 
nuclear mass. A multi-component fit using four elemental species suggests the 
same composition trends exhibited by the mean quantities, and also indicates that 
QGSJET and VENUS are the preferred hadronic interaction models. We find that 
an initially mixed composition turns lighter between 1 and 3 PeV, and then becomes 
heavier with increasing energy above 3 PeV. 

Key words: PACS 95.85.R. Cosmic rays. Knee, Energy spectrum. Composition, 
Cherenkov. 
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1 Introduction 



The all-particle energy spectrum of cosmic rays can be described by a steeply 
falling power law over many decades of energy. The smoothness of this drop in 
intensity with energy is broken by a change in index of the power law just above 
IQ-'^^eV. While the origin of this break (referred to as the "knee") is not yet 
fully understood, the prevailing theoretical models describe the knee as a result 
of the energy limit for particle acceleration in supernova shocks [1,2]. Further, 
these models predict that the composition of the primary cosmic rays should 
change from proton (or "light" -nuclei) dominated to iron (or "heavy" -nuclei) 
dominated as energy increases through the region of the knee. Measuring a 
composition trend from light to heavy would lend support to the supernova 
shock acceleration picture. 

Determining the composition of cosmic rays with energies greater than 10^^ eV 
is a notoriously difficult problem. To detect primary cosmic rays above this 
energy directly by satellite experiments requires an unacceptable launch pay- 
load volume. Similarly, stratospheric balloon-borne experiments are limited by 
volume and flight time in their collection of primary particles. Thus, to investi- 
gate the composition of cosmic rays at the knee, we must rely on ground-based 
detection of air showers generated by the primary cosmic rays. 

The Cherenkov light emission from the charged particle component of an 
air shower provides an integrated measurement of the longitudinal develop- 
ment [3,4] . One approach is to sample the Cherenkov lateral distribution, the 
photon density as a function of distance from the air shower core. The Cher- 
enkov intensity is proportional to the primary energy, while the slope of the 
lateral distribution is related to the depth of maximum shower development — 
and hence to the mass of the primary cosmic ray nucleus. Therefore measuring 
a large number of Cherenkov lateral distributions can provide information on 
how the composition changes with energy [5]. Previous attempts to exploit 
this fact at the knee include [6-8]. 

Optical photons suffer little absorption as they travel through the atmosphere. 
This means that the Cherenkov lateral distribution is much broader than that 
of charged particles. Additionally their numerical density is much higher. Thus 
it is possible to make high signal-to-noise measurements of Cherenkov lateral 
distributions using an array of detectors with smaller area and wider spacing 
than would be required for equivalent measurements of charged particles. 

To obtain high quality Cherenkov lateral distribution data the Broad Lat- 
eral Non-imaging Cherenkov Array (BLANCA) was built at the Chicago Air 
Shower Array (CASA) installation in Dugway, Utah. Using CASA as the cos- 
mic ray trigger, BLANCA operated on clear, moonless nights in 1997 and 
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Fig. 1. The C ASA-BLANC A arrays. The CASA scintiUator stations were placed 
at 15 m intervals. The 144 BLANCA air Cherenkov detectors filled the array with 
characteristic separations of 35-40 m. 

1998. In the following analysis we use CASA to find the shower core position 
and arrival direction and BLANCA to make a precision measurement of the 
Cherenkov lateral distribution. This paper details the results obtained through 
these measurements on the energy spectrum and composition of cosmic rays 
in the energy range between 3 x 10^"^ and 3 x 10^^ eV. 



2 The CASA-BLANCA Instrument 



The CASA-BLANCA instrument was located at the Dugway Proving Ground 
near Salt Lake City, Utah, USA, under a mean atmospheric overburden of 
870gcm-2. During BLANCA runs, CASA [9] consisted of 957 scintillation 
counters which detected the charged particles in an air shower. The surface 
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array covered approximately 0.2 km^. BLANCA [10] consisted of 144 angle- 
integrating detectors which recorded the lateral distribution of air shower 
Cherenkov light. The BLANCA detectors were not uniformly spaced but had 
an average separation of 35-40 m. MIA, an array of buried muon detectors at 
the same site, was not used in this analysis. Figure 1 shows the site plan. 

Each BLANCA detector contained a large Winston cone [11] which concen- 
trated the light striking an 880 cm^ entrance aperture onto a photomultiplier 
tube. The concentrator had a nominal half-angle of 12.5° and truncated length 
of 60 cm. The geometrical concentration ratio of this design was 19, but losses 
in the system reduced the effective concentration ratio to 15. Lab measure- 
ments and simulations show that the effective half-angle was actually ~ 10° 
because of a 6 mm gap between the photomultiplier and the cone. The Winston 
cones were aligned vertically with ~ 0.5° accuracy. A two-output preamplifier 
increased the dynamic range of the detector. A typical BLANCA unit had a 
detection threshold of approximately one blue photon per cm^. 

The Cherenkov array did not have a separate trigger system, relying instead 
on triggers from CAS A. Cherenkov data were recorded for all stations which 
exceeded a fixed threshold in coincidence with a surface array trigger. For 
showers in the BLANCA field of view, the CASA trigger threshold imposes a 
minimum energy of ~ 100 TeV on the Cherenkov array. However, to eliminate 
composition bias inherent in the CASA trigger threshold, we use only showers 
with an energy of at least 200 TeV as determined by BLANCA. 



3 Data Collection and Calibration 

CASA-BLANCA operated on 90 moonless nights between January, 1997 and 
May, 1998. After removing periods of hazy or cloudy weather, approximately 
460 hours of Cherenkov observations remain. Data cuts require events to 
have at least five good Cherenkov measurements from BLANCA and a recon- 
structed primary direction within 9° of zenith. Events are also cut if the core 
location reconstructed by CASA is outside the array or within 30 m of the edge, 
because core location uncertainty increases towards the edge. The geometrical 
and temporal cuts result in an exposure to cosmic rays of 1.83 x lO^^m^srs. 

For each night of data the BLANCA detectors are intercalibrated using the 
cosmic ray data itself to find their relative sensitivity to Cherenkov light. The 
intercalibration method depends on the circular symmetry of the Cherenkov 
light pool about the shower axis. Two detectors equidistant from a shower 
core receive, on average, equal Cherenkov photon densities, and any discrep- 
ancies can be attributed to different detector sensitivities. By averaging over 
suitable events in a run, we find the sensitivity ratio of each possible pair of 
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BLANCA detectors. A maximum-likelihood method is then used to determine 
the set of 143 relative gains that best reproduces the pairwise ratios [12]. The 
reliability of this method was verified with an in situ intercalibration using a 
portable stable blue LED flasher (A ~ 430 nm). The relative detector sensi- 
tivities are distributed log-normally with typical RMS of 0.4 in natural log. 
The relative sensitivies are stable over time except for occasional changes in 
detector hardware. 

The BLANCA absolute calibration also used the blue LED system as a refer- 
ence source. The LED flasher was calibrated by using it to produce single pho- 
toelectrons in a photomultiplier with good charge resolution. Two BLANCA 
detectors were then calibrated in a dark box using the reference LED. The 
absolute calibration has a 20% systematic uncertainty, resulting mainly from 
the uncertain quantum efficiency of the reference photomultiplier. This ab- 
solute calibration error produces a similar uncertainty in the overall energy 
scale. The spectral response of the BLANCA photomultipliers was measured 
by the manufacturer for four representative tubes. As the relative response 
was similar for these four detectors, it was also assumed to be the same for all 
BLANCA detectors. 



4 Air Shower Simulations 

The CASA-BLANCA analysis compares the Cherenkov measurements with 
air showers simulated by the CORSIKA Monte Carlo version 5.621 [13]. We 
used the EGS4 and GHEISHA codes for the electromagnetic and low-energy 
nuclear interactions. Nucleus-nucleus interactions at air shower energies are 
well beyond the reach of accelerator experiments. Therefore we are forced to 
rely on hadronic interaction models which attempt to extrapolate from the 
available data using different mixtures of theory and phenomenology. Several 
groups produce such models — in this paper we have used four: QGSJET [14], 
VENUS [15], SIBYLL [16], and HDPM [13]. BLANCA data are interpreted 
according to each model, indicating the systematic errors that depend on the 
choice of interaction model. 

Simulated showers were produced for proton, helium, nitrogen, and iron pri- 
maries in equal numbers. Nitrogen was chosen to represent the entire CNO 
group. The shower libraries consist of 10,000 showers per primary species, per 
hadronic model. Primary energies are distributed uniformly in log(£') between 
10^^ and 10^^'^ eV. Shower directions are uniform in solid angle to a maximum 
zenith angle of 12°. To produce such large simulated air shower libraries, it 
was necessary to employ the CORSIKA thinning option, which tracks only 
a representative sample of shower particles below a threshold energy. In all 
simulations, this threshold was 10~^ times the primary energy. Studies of sim- 
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ulated 1 PeV showers showed that at this thinning level, the distributions of 
gross properties such as the depth of shower maximum and Cherenkov slope 
and intensity were indistinguishable from those of unthinned showers [17]. 

Shower development and Cherenkov emission were simulated in Monte Carlo 
assuming the U. S. standard atmosphere. The CORSIKA program was mod- 
ified to include a complete model of atmospheric scattering, both Rayleigh 
(molecular) and Mie (aerosol). Although the modified CORSIKA tracks scat- 
tered photons, few reach the ground within the BLANCA field of view, effec- 
tively making scattering an absorption process for Cherenkov photons. The 
scattering losses were generally similar in magnitude to the expected measure- 
ment errors and were correlated with the depth of shower maximum. There- 
fore, atmospheric scattering must not be ignored in analyzing the BLANCA 
data. On average, scattering reduces the Cherenkov intensity by ~ 20% and 
increases the inner slope by ~ 7%. Scattering effects are smallest for late 
developing showers. Possible molecular absorption by oxygen and ozone was 
determined to be small compared with the systematic error in the absolute 
detector calibration and was consequently ignored. 

The CORSIKA air showers were processed by a full BLANCA detector sim- 
ulation which includes the measured wavelength dependence and angular re- 
sponse of the BLANCA detectors. Other simulated effects include detector 
alignment; unequal detector gains and saturation levels; night sky background 
light; photomultiplier response; and errors in the CAS A core location and 
shower direction. The detector simulation produces "fake data" which is cal- 
ibrated and fit like the real data. We have used this fake data to find the 
optimum transfer functions for converting Cherenkov measurements to air 
shower and primary cosmic ray parameters throughout this paper. 



5 The Cosmic Ray Energy Spectrum 

For each air shower event, raw BLANCA data are converted to photon densi- 
ties, producing a Cherenkov lateral distribution. We fit this lateral distribution 
with an empirically motivated function which matches both the real and sim- 
ulated data. The function is exponential in the range 30 m-120 m from the 
shower core and a power law from 120m-350m. It has three parameters: a 
normalization C120, the exponential "inner slope" s, and the power law index 

\ C'120 e"(i20 m-r) 30m<r<120m 
C{r) = \ - (1) 

[ C120 (r/120m)-^, 120m < r < 350m 
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The energy of each air shower is derived using only the Cuo and s parameters 
of the Cherenkov lateral distribution fit. The outer slope /? is not used, both 
because it correlates strongly with s and because it is subject to larger mea- 
surement errors. The Monte Carlo fake data hbraries (including detector sim- 
ulation) are used to determine the relationship between measured quantities 
and energy. The energy depends primarily on C120, the Cherenkov intensity 
120 m from the core. We fit the logarithm of the energy as a quadratic function 
of logCi2o (the curvature is small, typically 0.005 decades"^). In all hadronic 
models, C120 grows approximately as because the fraction of primary 

energy directed into the electromagnetic component of the cascade increases 
with energy. 

According to the Monte Carlo, the quadratic function used to estimate the 
energy works well for most showers but has a bias such that energies are 
slightly underestimated for showers with unusually large or small depths of 
maximum development. Therefore, a small correction is applied to the energy 
estimate. The correction depends on the Cherenkov slope s and on the shower 
zenith angle. The magnitude of this energy correction is less than 10% for 85% 
of BLANCA showers. 

Energies derived from the data in this manner have an error distribution which 
depends on the primary mass and energy. In general, the random errors on 
reconstructing a single shower's energy are comparable to the systematic un- 
certainty due to the unknown composition. Assuming a mixed cosmic ray 
composition, the BLANCA energy resolution for a single air shower is approx- 
imately 12% for a 200 TeV shower, falling to 8% for energies above 5 PeV. 

The differential all-particle cosmic ray flux measured by CASA-BLANCA is 
shown in Figure 2, scaled up by a factor of (ii^/1 GeV)^'^^ to emphasize the 
structure. Table A.l lists the values of the observed spectrum. The energy 
spectrum is compared with that reported by several other groups. Although 
the CASA-BLANCA, DICE, and CASA-MIA experiments shared some instru- 
mentation, their data sets and energy analysis methods are entirely indepen- 
dent. These three experiments show very good agreement in their spectrum 
determination. Most other results are consistent with CASA-BLANCA, given 
the 20% or larger energy systematic error typical of air shower measurements. 

The spectrum shown in Figure 2 uses event energies derived from the Cheren- 
kov predictions of CORSIKA with the QGSJET hadronic interaction model. 
The data can also be interpreted using the other available interaction models. 
The alternate energy estimates lead to spectra with no qualitatively different 
features. Instead, they amount only to a shift in energy scale of order 10%, 
less than the BLANCA instrumental energy scale uncertainty. HDPM and 
VENUS predict less Cherenkov light and hence assign higher energies than 
QGSJET does, while the SIBYLL simulations lead to lower assigned energies. 
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Fig. 2. The differential all-particle cosmic ray flux measured by C ASA-BLANC A 
(QGSJET energy). Each data point represents the differential flux scaled up by a 
factor of GeV)^'^^; the error bars represent statistical (Poisson) errors only. 
The diagonal bar shows the effect of a possible systematic energy shift of ±18%. 
Also plotted are the cosmic ray fluxes reported by DICE [18], CASA-MIA [19], 
Akeno [20] and Tibet AS7 [21]. 



Figure 2 contains the knee region of the spectrum, near 3 PeV. The CASA- 
BLANCA spectrum exhibits a smooth change rather than a sharp break here. 
However, measurements over a wider energy range show that the form of the 
cosmic ray spectrum is a power law well above and below 3 PeV. Historically, 
many groups have found the knee to be quite sharp. In the spirit of the usual 
discussion of the knee, we have fit several similar functions to the data (Fig- 
ure 3, top). The fits find simultaneously the position (energy) of the knee and 
the power law indices above and below the knee. A log-likelihood fit is per- 
formed in order to account for the Poisson statistics of discrete events in a 
binned energy distribution. 
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Fig. 3. Top: Several trial functions used to fit the cosmic ray flux. Bottom: The 
log-likelihood for the smooth function at a range of widths w. The CASA-BLANCA 
data fit best to a knee with w = 0.25, or a full width of one-half decade. 

Of the trial functions, the smoothly changing power law fits the BLANCA 
data best: 



J{E) = Jk 



E 
Eh 



1 + 



E 



{l3—a)w 



(2) 



Ek is the energy at the center of the transition, i.e. the knee energy. For 
E <^ Ek, the function is a power law with index a, while the spectral index 
becomes (3 lor E ^ E^. Parameter sets the normahzation at the knee. The 
fifth parameter, is the half-width in decades of the transition region. The 
lower panel of Figure 3 shows how the log-likelihood depends on the choice 
of w. The data favor a knee one-half decade wide {w = 0.25). The best-fit 
knee energy is 2.0loi PeV, with power law indices of a = —2.72 ± 0.02 and 
P — —2.95 ± 0.02. Figure 4 shows the energy spectrum near the knee using 
fine bins 0.04 decades wide as well as the best fit curve. 
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Fig. 4. The CASA-BLANCA energy spectrum shown with more detail in the knee 
region. 

6 Depth of Shower Maximum 

The depth of shower maximum (X^ax) is an important characteristic of air 
shower development which, for given energy, is related to the mass of the 
primary particle. However, like all air shower parameters, the relationship 
depends on the choice of uncertain high energy hadronic interaction model 
parameters. Imaging experiments such as Fly's Eye [22], and to a lesser extent 
DICE [18], measure Xmax rather directly. The importance of Xmax for various 
other types of ground-based experiment has long been known [23]; measured 
parameters can often be translated into the depth of shower maximum in a way 
which is rather independent of the hadronic model. Therefore X^ax provides 
a useful middle ground on which experiments may publish and compare their 
results. 

The mean Xmax for a given primary type grows logarithmically with energy at 
an approximate elongation rate of 80gcm~^ per decade, although this value 
depends on the hadronic model used. The expected Xmax is similar for two 
primaries of different mass if they have equal energy per nucleon. 

To determine the optimum transfer function for converting Cherenkov lateral 
distributions into Xmax, we study the same set of simulated showers used 
to derive the primary energy function. The fake data libraries use the four 
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Fig. 5. The mean depth of shower maximum {Xmax) measured by C ASA-BLANC A 
as a function of energy. The hnes represent the values for pure proton or pure 
iron samples predicted by each of the four hadronic interaction models. The thick 
error bars represent the statistical uncertainty, while the thin errors also include the 
systematic uncertainties discussed in the text; the former are important only above 
lOPeV. 

primary types and four hadronic interaction models processed through the 
BLANCA detector simulation. The simulated Cherenkov lateral distributions 
are fit to the function in Equation 1, which is exponential with a slope s from 
30 m to 120 m from the core. The combined shower and detector simulations 
show that this inner slope is linearly related to X^ax except for the deepest 
developing showers; an additional small quadratic term is required for slopes 
exceeding = 0.018 m~^. 

The mean X^ax is shown in Figure 5 as a function of energy. Both quantities 
are derived from the CASA-BLANCA data using the CORSIKA/QGSJET 
Monte Carlo results. We indeed find that the results are very similar if any 
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Fig. 6. Systematic uncertainties on CASA-BLANCA X^aax estimates. At low ener- 
gies, the dominant error comes from the slope-to- X^ax conversion function; at high 
energies, the uncertain photomultiplier linearity is more important. For compari- 
son, the statistical error on the CASA-BLANCA mean X^aax measurements is also 
shown (dashed line), assuming a bin width of 0.1 decades. 



of the other hadronic models are used instead. Numerical results are given in 
Table A. 2. Figure 5 also shows the mean X^ax expected for pure samples of 
proton primaries and iron primaries. SIBYLL generally predicts deeper shower 
maximum than other models, while HDPM exhibits a steeper elongation rate 
than the others (~ 90gcm~^ per decade compared with 70-75 gcm~^ typical 
of the other models). The BLANCA results are clearly consistent with a mixed 
composition throughout the energy range, regardless of the preferred hadronic 
model. The data suggest that the composition becomes lighter approaching 
the knee and then becomes heavier at higher energies. 

The shower depth estimated from Cherenkov observations is subject to a num- 
ber of small systematic uncertainties. Random core errors lead to a systematic 
flattening of the Cherenkov inner slope, but the effect is only a very small bias 
toward deeper X^ax- Photomultiplier saturation poses a potential problem at 
high energies. However, the nonlinearity has been characterized in laboratory 
studies of fourteen BLANCA detectors and the data corrected for its effects. 
The uncertainty on this correction leads to a systematic error in X^ax which is 
only lOgcm"^ at the highest energies. A third error dominates at energies be- 
low 1 PeV. This error arises from the limitations of the function which converts 
Cherenkov slope to an estimate of X^ax- The function tends to overestimate 
the depth of showers at the extreme ends of the BLANCA energy range. It 
is difficult to overcome this weakness without introducing at the same time 
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a much larger bias which depends on X-max itself. Instead, we take the error 
found in Monte Carlo studies as a systematic error on the measured X^ax- 
The independent errors are added in quadrature to find the total systematic 
error (Figure 6). Systematic errors are important only below lOPeV. At high 
energy, statistical errors are the more serious limitation on measuring mean 



7 Mean Nucleeir Mass 



The mean depth of shower maximum results presented above are essentially 
independent of a particular hadronic interaction model. They are therefore 
useful for comparison with other experiments and re-interpretation on the 
basis of future hadronic interaction models. However, a quantity which is of 
much more direct astrophysical interest is the mean nuclear mass of the pri- 
mary cosmic rays. 

We choose to derive mean primary mass directly from the Cherenkov lateral 
distribution slope s. It would make little difference if we were to do so via 
the Xmax values discussed in the previous section. The important point is 
that while the transfer function from s to X^ax is rather independent of the 
hadronic interaction model, any interpretation in terms of absolute nuclear 
mass is not. This is clear from the disparity among the models of the mean 
proton and iron X^ax values shown in Figure 5. 

At fixed primary energy, s and X^ax both depend linearly on the logarithm of 
nuclear mass A, as do most composition-sensitive air shower parameters. Fol- 
lowing previous authors, we choose to work with the natural log, ln(A). Unlike 
Xmax=fi{s), the transfer function ln(/l)=/2(s) depends on energy. Therefore 
we divide the Monte Carlo fake data into six bands of Cuo and perform a 
linear fit to \ii{A) versus s in each. To interpret each real event the slope and 
intercept of the transfer function are interpolated between the appropriate 
bracketing bands. 

This method has little systematic bias apart from the differences between 
hadronic interaction models. The mean reconstructed \n{A) for a pure sample 
of each simulated primary species is accurate to to ±20% over the BLANCA 
energy range. On the other hand, the random error is large on any single 
measurement of ln(74). Shower-by-shower estimates of primary mass cannot 
be made with any accuracy — the fluctuations inherent in the air shower 
process preclude them. Nevertheless, the mean value of In(^) provides a useful 
indicator of the cosmic ray mass composition. 

The mean In (A) is shown as a function of primary energy in Figure 7. The four 
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Fig. 7. The mean logarithmic mass (ln(A)) measured by C ASA-BLANC A as a 
function of energy. The four sets of symbols show the BLANCA data interpreted 
using CORSIKA coupled with the indicated hadronic interaction model. Error bars 
are statistical only and shown only on the QGSJET results; errors on the other 
points are of similar size. Assuming a hadronic model, the systematic uncertainties 
in the ln(A) estimate are typically 0.2. The dominant systematic effect is clearly 
the difference between the hadronic models. 

sets of symbols show the BLANCA data interpreted using CORSIKA coupled 
with each of the hadronic codes. Numerical results are given in Table A. 3. The 
dependence on interaction model is clear: the models set the overall mass scale 
differently, but they indicate the same mass variation with energy. The mean 
mass becomes lighter with increasing energy through the knee, then becomes 
heavier above ~ 3 PeV. The ln(y4) plot exhibits the same trends seen in the 
^max results of the previous section. 
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8 A Multi-species Fit to the Cherenkov Data 



The techniques described above involve estimating X^ax or \ii{A) for each 
shower and then taking the average over all showers in a given energy range. By 
considering only the average value we lose much of the available information. 
Instead the measured distribution of a composition sensitive parameter can 
be compared with those predicted for a number of simulated primary species, 
providing a more powerful technique to study cosmic ray composition. 

Comparing measured and simulated distributions requires high statistics sam- 
ples for both the real and Monte Carlo data sets. We separate the data into five 
logarithmic energy bins between 10^"^'^ and 10^^'^ eV. This bin choice is a com- 
promise between the need to include many showers in each range and the wish 
to examine trends on as fine an energy scale as possible. Within each range, 
we find the distribution of the Cherenkov inner slope (s) for the real data and 
for pure samples of each species in the Monte Carlo library (protons. He, N, 
and Fe). The simulated distributions of s (with detector effects) are smoothed 
by a multiquadratic smoothing algorithm [24] . To preserve information about 
their limited statistics, the data distributions are not smoothed. 

The multi-species fit in each energy range finds the linear combination of the 
four simulated distributions which reproduces the data distribution best. Since 
each primary species has a characteristic shape of its s distribution, this fit 
uses more information than simply the mean or even the width of s. We do not 
a priori require the fractional contribution of each primary type to lie in the 
physical range of 0-100%, nor is the sum of the fractions constrained to equal 
1.0. In practice, however, the sum is always in the range 100 ± 0.5%. The 
fits use a MlNUlT-based log-likelihood maximization procedure [25], which 
accounts properly for the Poisson probability distribution of data events in 
bins with low statistics. 

As an example. Figure 8 shows the multi-species fit in the energy range 10^^'^- 
-[^q15.3 gY r^Yie lower right panel (#4) displays the full fit using all four available 
primary types. The other panels show the best fits that can be made when 
hehum (#2), nitrogen (#3), or both (#1) are omitted. The Monte Carlo pre- 
dictions cannot match the data using protons and iron alone; the intermediate 
mass nitrogen species is also required. Panel ^2 shows that the best fit of p, N, 
and Fe to the data is very close but fails to match the shape at the peak and 
in the long tail to high values of s. The data strongly suggest that at least the 
four primary types considered here contribute to the cosmic rays just below 
the knee. 

The fits demonstrated in Figure 8 (lower right) are performed for all five energy 
ranges and using the predictions of all four high energy hadronic interaction 
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Fig. 8. A demonstration of the multi-species fitting procedure, in this case the energy 
range 

1014.9 _ 1015-3 gv. The sohd histo gram gives the distribution of Cherenkov 
slope s for all CASA-BLANCA data in that energy range and is the same in all four 
panels. The solid curves show the combinations of Monte Carlo (QGSJET) showers 
which best reproduce the measured distribution. Panels 1-4 correspond to different 
combinations of primary species. Clearly a nitrogen component is needed to match 
the data, but a helium component is also important. 



models. Gauging the goodness of fit presents a problem. As stated above the 
multi species fit uses a log-likelihood maximization procedure; although this 
is an appropriate method for extracting abundance fractions from the binned 
data, the actual value of the likelihood is not useful [25]. Therefore we calculate 
and use the familiar quantity as a rough guide to the fit quality. In the lowest 
energy range, the number of data events is so large that no combination of the 
four cosmic ray species can reproduce the data adequately. This is probably 
the result of limitations of the shower and detector simulations, although it 
could also be due to the limited number of species considered. Conversely, the 
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Fig. 9. Results of the multi-species fit to the C ASA-BLANC A data. The line graphs 
indicate the mixture of proton, helium, nitrogen, and iron primaries which best 
reproduces the Cherenkov slope distributions of the data. The lines show cumulative 
fractions (i.e. the lowest line gives the proton fraction, while the next line gives 
the combined proton and helium fraction). QGSJET (left) and VENUS {right) 
are shown. The other two interaction models show similar trends but with heavier 
overall composition. The stars (*) at 10^^ eV show the results of direct measurements 
from the JACEE balloon-borne emulsion experiment [26,27]. JACEE Ne-Si data 
have been divided evenly into the N and Fe groups for comparison with BLANCA. 



high energy ranges have too few events to constrain the abundances well. The 
results for all models are presented in Table A. 4. 

The results of the multi-species fit to the BLANCA Cherenkov slope data are 
shown in Figure 9 for the QGSJET and VENUS models. The SIBYLL and 
HDPM models show similar trends but a heavier overall composition. These 
latter two models also give unphysical negative helium abundances in at least 
one energy bin, and systematically poorer fits. At 100 TeV, data from the 
JACEE balloon direct measurements [26,27] are shown for comparison. The 
direct composition at 100 TeV agrees well with the BLANCA data at 400 TeV. 
The results of the multi-species fit also agree with the mean X^ax and mean 
ln(A) derived in the previous two sections. All three ways of interpreting the 
data indicate that the cosmic ray composition is lighter near 3 PeV than it is 
at either 300 TeV or 30 PeV. 
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9 Conclusions 



The CASA-BLANCA experiment has studied cosmic rays in the energy range 
0.3-30 PeV. The primary energy and mass are found by measuring the Cher- 
enkov lateral distribution for each air shower. In an effort to understand how 
results depend on the unknown physics of high energy nuclear interactions, 
we have interpreted the data using the CORSIKA air shower Monte Carlo 
program with four different hadronic interaction models: QGSJET, VENUS, 
SIBYLL, and HDPM. 

The BLANCA energy spectrum agrees well with previous measurements and 
exhibits a smooth knee near 2-3 PeV in primary energy. The model dependence 
of the energy scale is less than the absolute calibration uncertainty. 

We find the transformation from measured Cherenkov lateral distribution 
slope to the depth of shower maximum Xmax to be essentially model inde- 
pendent. In Figure 10 our results are compared to previous experiments over 
a wide energy range. The BLANCA data are well within the physically rea- 
sonable range bounded by the pure proton and iron curves; furthermore they 
are consistent at low energy with those expected from direct measurements 
and at high energy with the Fly's Eye result [28] . 

We have also interpreted our data as a mean nuclear mass. This is essentially 
equivalent to the X^ax analysis but is a quantity of more direct astrophysical 
interest. 

It has been a long held goal in the air shower field to choose an adequate 
hadronic interaction model and determine the nuclear composition of the pri- 
mary cosmic rays simultaneously. With the advent of the powerful simulation 
tool provided by the CORSIKA group, and high collecting power arrays such 
as CASA-BLANCA, it seems that this ambition may be becoming a reality. A 
multi component fit of the type described in Section 8 is a much more efficient 
use of the available data than simply considering the mean value and spread of 
a quantity, and the experimental statistics are starting to justify this approach. 
The agreement between data and simulation in Figure 8 is impressive. 

On the basis of our data we favor the QGSJET and VENUS models and reject 
SIBYLL and HDPM. At the same time, both of the ways in which we have an- 
alyzed our data indicate that the cosmic ray composition is lighter near 3 PeV 
than it is at either 300 TeV or 30 PeV. The trend towards heavier primary 
mass above 3 PeV agrees with the canonical model of Galactic production and 
a rigidity-dependent time for escape, and is not consistent with acceleration 
at sites such as AGN, which require a pure proton composition well above the 
knee [30]. 
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Fig. 10. The CASA-BLANCA measurement of X^ax compared with other results. 
All experiments operating near the knee use atmospheric Cherenkov light, including 
DICE [18], the AIROBICC array of HEGRA [7], and the VULCAN array at the 
South Pole [8]. The high energy measurements (> 10^^ eV) by Fly's Eye use the 
atmospheric fluorescence technique [28]. The "direct" point estimates the mean 
Xmax that would be expected on the basis of direct balloon measurements [18]. The 
Monte Carlo hues use CORSIKA with QGSJET [29]. 

The trend shown in our data to a lighter composition approaching the knee 
is puzzling but not without theoretical precedent. Swordy, arguing that there 
must be a minimum path length in the Galaxy even for the highest energy 
cosmic rays, predicts a light composition at the knee [31]. 
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A Data tables 



Energy range Differential flux J{E) in each bin 
logio(-E/leV) (m-^sr-is-^GeV-i) 



14.3 - 


14.4 


(12.12 ±0.04) X 10-1^ 


14.4 - 


14.5 


(6.68 ± 0.02) 


X 10-11 


14.5 - 


14.6 


(3.58 ± 0.02) 


X 10-11 


14.6 - 


14.7 


(1.91 ±0.01) 


X 10-11 


14.7 - 


14.8 


(1.03 ±0.01) 


X 10-11 


14.8 - 


14.9 


(5.42 ± 0.04) 


X 10-12 


14.9 - 


15.0 


(2.90 ± 0.03) 


X 10-12 


15.0 - 


15.1 


(1.57 ±0.02) 


X 10-12 


15.1 - 


15.2 


(8.18 ±0.12) 


X 10-13 


15.2 - 


15.3 


(4.36 ± 0.08) 


X 10-13 


15.3 - 


15.4 


(2.21 ± 0.05) 


X 10-13 


15.4 - 


15.5 


(1.22 ±0.03) 


X 10-13 


15.5 - 


15.6 


(6.2 ±0.2) X 10-1^ 


15.6 - 


15.7 


(2.9 ±0.1) X lO-i"^ 


15.7- 


15.8 


(1.5 ±0.1) X 10"^^ 


15.8 - 


15.9 


(7.7 ±0.5) X 10-15 


15.9 - 


16.1 


(2.9 ± 0.2) X 10-15 


16.1 - 


16.3 


(8.1 ± 0.8) X 10-16 


16.3 - 


16.5 


(2.1 ±0.3) X 10-16 


16.5 - 


16.7 


(3.1 ± 1.0) X 10-1^ 


16.7 - 


16.9 


(2.3 ± 0.7) X 10-1^ 



Table A.l 

The primary cosmic ray energy spectrum measured by CASA-BLANCA. Bin widths 
rise with increasing energy so that Emax/Emin = lO'^'i at lower energies, while 
Emax/Emin = 10° 2 for the five highest bins. Errors represent only the Poisson 
uncertainty in each bin. There is an additional instrumental systematic uncertainty 
of 18%. These results use the QGSJET-derived energy transfer function. 
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Energy range Mean s 
logio(£^/eV) (10-3 m-i) 



i^max) ± Stat. ± sys. 
(gcm-2) 



14.3 - 


14.4 


11.4 ±0.0 


458 ± 0.3 ± 12 


14.4 - 


14.5 


11.9 ±0.0 


469 ± 0.3 ± 10 


14.5 - 


14.6 


12.5 ±0.0 


479 ± 0.4 ± 9 


14.6 - 


14.7 


12.9 ±0.0 


488 ± 0.4 ± 9 


14.7 - 


14.8 


13.4 ± 0.0 


497 ± 0.5 ± 9 


14.8 - 


14.9 


13.7 ±0.0 


504 ± 0.6 ± 9 


14.9 - 


15.0 


14.1 ±0.0 


511 ±0.8 ±8 


15.0 - 


15.1 


14.7 ±0.0 


523 ± 1 ± 7 


15.1 - 


15.2 


15.2 ±0.1 


532 ± 1 ± 7 


15.2 - 


15.3 


15.6 ±0.1 


542 ± 1 ± 7 


15.3 ~ 


15.4 


16.3 ±0.1 


555 ± 2 ± 7 


15.4 - 


15.5 


16.8 ±0.1 


565 ± 2 ± 7 


15.5 - 


15.6 


17.2 ±0.1 


573 ± 3 ± 7 


15.6 - 


15.7 


17.3 ±0.2 


576 ± 4 ± 8 


15.7- 


15.8 


17.4 ±0.2 


578 ± 4 ± 8 


1 c; Q 

10. o — 


10. y 


1 7 tc -I- n 
1 ( .0 ± u.z 


/ y ± ± o 


15.9 - 


16.0 


17.9 ±0.3 


588 ± 6 ± 9 


16.0 - 


16.1 


17.0 ±0.4 


570 ± 8 ± 9 


16.1 - 


16.2 


17.5 ±0.4 


579 ± 10 ± 9 


16.2 - 


16.3 


17.4 ±0.6 


576 ± 12 ± 10 


16.3 - 


16.4 


18.1 ±0.5 


589 ± 11 ± 11 


16.4 - 


16.5 


18.5 ±0.8 


600 ± 19 ± 11 


16.5 - 


16.6 


17.2 ± 1.5 


570 ± 31 ± 12 



Table A. 2 

The mean Cherenkov inner slope (s) measured by CASA-BLANCA and the corre- 
sponding mean Xmax- The first column of errors given for each quantity is statistical; 
the standard deviation divided by ^/N. For X,nax a systematic error is given which 
is due to a combination of effects (see Section 6). These results use the QGSJET- 
derived X^ax transfer function. 
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Energy range Mean In(^) 

logio(^/leV) QGSJET VENUS SIBYLL HDPM 



14.3 - 


14.4 


1.71 ± 0.01 


1.58 


2.27 


1.99 


14.4 - 


14.5 


1.67 ± 0.01 


1.63 


2.24 


2.06 


14.5 - 


14.6 


1.63 lb 0.01 


1.69 


2.28 


2.14 


14.6 - 


14.7 


1.61 ± 0.01 


1.71 


2.29 


2.18 


14.7 - 


14.8 


1.60 ± 0.02 


1.74 


2.35 


2.20 


14.8 - 


14.9 


1.69 ± 0.02 


1.80 


2.47 


2.26 


14.9 - 


15.0 


1.65 ± 0.03 


1.85 


2.45 


2.32 


15.0 - 


15.1 


1.49 ± 0.03 


1.78 


2.34 


2.32 


15.1 - 


15.2 


1.39 ± 0.04 


1.69 


2.29 


2.18 


15.2 - 


15.3 


1.28 ± 0.05 


1.61 


2.16 


2.10 


15.3 - 


15.4 


1.07 lb 0.06 


1.36 


1.99 


2.01 


15.4 - 


15.5 


jA rt rt 1 rt /A^^ 

0.90 lb 0.07 


1.28 


1.78 


1.83 


15.5 - 


15.6 


f\ r\ 1 /A /Art 

0.86 lb 0.09 


1.06 


1.68 


1.59 


15.6 - 


15.7 


rt rt K 1 rt -1 rt 

0.95 lb 0.12 


1.25 


1.81 


1.85 


15.7 - 


15.8 


1.1 lb 0.1 


1.2 


2.0 


1.7 


15.8 - 


15.9 


1.2 lb 0.2 


1.6 


2.0 


2.0 


15.9 - 


16.0 


1.0 lb 0.2 


1.3 


2.1 


2.1 


16.0 - 


16.1 


2.0 zb 0.3 


2.3 


2.9 


2.5 


16.1 - 


16.2 


1.8 ±0.4 


2.1 


2.8 


2.7 


16.2 - 


16.3 


2.1 ±0.5 


2.3 


2.6 


2.9 


16.3 - 


16.4 


1.7 ±0.5 


1.9 


3.1 


2.4 


16.4 - 


16.5 


1.4 ±0.8 


2.2 


2.9 


3.2 


16.5 - 


16.6 


2.9 ± 1.4 


3.2 


2.8 


2.4 



Table A.3 

The mean In(^) measured by CASA-BLANCA. The four columns show the data 
interpreted according to each hadronic interaction model. The variation provides 
some insight into the systematic errors. Statistical errors shown on the QGSJET 
points are similar for all four models. 
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Energy Range Abundance (%) of 



logio(-E/leV) p He N Fe Fit 









QGSJET 








14.5 - 


14.9 


21.8 ±0.4 


40.1 ±0.7 


23.4 ± 0.7 


14.6 ± 0.3 


38.1 


14.9 - 


15.3 


42 ± 1 


21 ±2 


19 ± 1 


18 ± 1 


4.5 


15.3 - 


15.7 


51 ±3 


33 ±4 


3±3 


13 ± 1 


1.9 


15.7- 


16.1 


53 ±8 


14± 10 


23 ±6 


10 ±3 


0.7 


16.1 - 


16.5 


31± 12 


12± 18 


35± 17 


22 ±8 


1.9 








VENUS 








14.5 - 


14.9 


23.9 ±0.4 


27.6 ±0.7 


31.8 ±0.5 


16.7 ±0.3 


47.8 


14.9 - 


15.3 


29 ± 1 


29 ±2 


19 ± 1 


23 ± 1 


5.9 


15.3 - 


15.7 


46 ±2 


23 ±4 


15 ±3 


16 ± 1 


1.7 


15.7- 


16.1 


46 ±6 


6±9 


33 ±7 


15 ±3 


0.8 


16.1 - 


16.5 


16 ±9 


33± 14 


23± 13 


29 ±8 


1.8 








SIBYLL 








14.5 - 


14.9 


16.8 ±0.4 


20.9 ±0.7 


24.8 ±0.7 


37.4 ±0.3 


49.4 


14.9 - 


15.3 


35 ± 1 


-7±2 


39 ± 1 


33 ±1 


21.6 


15.3 - 


15.7 


37 ±3 


21 ±5 


9±4 


33 ±2 


4.8 


15.7- 


16.1 


31 ±6 


19 ±10 


9±8 


41 ±4 


1.2 


16.1 - 


16.5 


9±9 


31 ±17 


8 ±19 


53± 11 


2.6 








HDPM 








14.5 - 


14.9 


19.9 ±0.3 


17.9 ±0.5 


31.4 ±0.5 


30.7 ±0.2 


90.0 


14.9 - 


15.3 


19 ± 1 


23 ± 1 


24 ± 1 


34 ±1 


11.5 


15.3 - 


15.7 


32 ±2 


16 ±3 


26 ±2 


26 ±1 


3.2 


15.7- 


16.1 


37 ±4 


-3 ±6 


43 ±6 


23 ±3 


1.2 


16.1 - 


16.5 


12 ±6 


21 ±12 


18± 12 


49 ±8 


1.5 



Table A.4 

Results of the multi-species fits to the CASA-BLANCA data. Statistical errors on 
each fraction are strongly correlated. The errors increase with energy due to limited 
statistics. Unphysical negative abundances are a result of a poor hadronic model 
and/or inadequate statistics. 
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